/******************************************************************************
Safer in School? The Impact of Compulsory Schooling on Maltreatment and Associated Harms
by Adam A. Dzulkipli, Nicole Black, David W. Johnston, and Leonie Segal

Description: This program recreates all of the Appendix tables presented in the
supplemental Online Appendix available for the paper. Tables are exported as .txt 
files using the 'esttab' package in Stata. If a table has multiple panels, the 
panels are appended together using the 'append' option in 'esttab'. The table 
will then require further manual editing in a typesetting software to get the 
table to the exact format as it is presented in the paper. The 'esttab' command
can also be altered to export the tables as .tex files if so desired.

Note that directory paths have been removed to preserve privacy. Please insert 
the relevant directories in the global macros under "Setup" before running the
.do file. This code was written and implemented in Stata MP 19.5, with
a computational time of approximately ~2.5 hours.

Packages required include:
•	ftools
•	reghdfe
•	estout
•	rdrobust
•	pzms
•	rddisttestk (.ado file by Brigham Frandsen available from
	https://sites.google.com/view/brighamfrandsen/software?authuser=0)
•	rddensity
•	rdhonest

*****************************************************************************/

********************************************************************************
*	Setup
********************************************************************************
clear all
set more off
macro drop all
capture log close

graph set window fontface "Times New Roman"
set scheme s2color

global raw 
global tmp 
global clean 
global output 

********************************************************************************
*	Appendix Table A1 to Appendix Table A3
********************************************************************************
/*
  These tables were constructed manually in LaTeX. The grouping of ICD-10 diagnosis 
  codes into injuries, mental illness, and digestive system issues are partially based
  from the classification of ED visits in Gnanamanickam et al. (2022).
*/

********************************************************************************
*	Appendix Table A4: Composition of ED visits by diagnostic blocks and CPS contact
********************************************************************************

use "${tmp}ed_detailed_intermediate.dta", clear

* Merge to final analytical file and keep only the relevant sample of children
merge m:1 pslk3 using "${clean}analysis_final.dta"
keep if _merge == 3
drop _merge

/*
  Construct category for missing Major Diagnostic Block (MDB) code (i.e. 
  unstated diagnosis)
*/
replace e_mdb_gr = "Unstated" if e_mdb_gr == "" & rectype == 1

/*
  Group drug/alcohol abuse and mental disorders (code 1D) with other ED visits 
  relating to mental health (code 1)
*/
replace e_mdb_gr = "1D" if e_mdb_gr == "1"

* Construct table of absolute and relative frequencies of ED visits
eststo clear

forval y = 0/1{
	foreach x in edvisit16 edvisit_17to21{
		 eststo: estpost tab e_mdb_gr if `x' == 1 & cps_age16_binary == `y', sort
	}
}
		
* Export table using the 'esttab' package
esttab using "${output}Appendix_TableA4.txt", label replace b(%9.4f) se(%9.4f) noomit ///
		cells("b pct(fmt(1))") obslast ///
		collabels("Frequency" "Percent") ///
		coeflabel(1D "Mental health" 2A "Injuries" ///
		3A "Circulatory system illness" 3B "Respiratory system illness" ///
		3C "Digestive system illness" 3D "Urological system illness" ///
		3E "Neurological system illness" 3F "Illness of the eyes" ///
		3G "Illness of the ear/nose/throat" ///
		3H "Musculoskeletal/connective tissue system illness" ///
		3I "Illness of skin/subcutaneous tissue/breast" ///
		3J "Blood/immune system illness" 3K "Obstetric illness" ///
		3L "Gynaecological illness" 3LM "Reproductive system illness" ///
		3M "Male reproductive system illness" 3N "System infection/parasites" ///
		3O "Illness of other/unknown systems" 3P "Newborn/neonate illness" ///
		3Q "Hepatobiliary system illness" ///
		3R "Endocrine/nutritional/metabolic system illness" 3S "Allergy" ///
		5 "Social problem" 6 "Other") ///
		nobase nogaps nomtitle nonumber float varwidth(50) ///
		order(2A 1D 3C) drop(Total) ///
		title("Appendix Table A4")

********************************************************************************
*	Appendix Table A5: Descriptive statistics for sample born before reform cut-off
********************************************************************************

use "${clean}analysis_final", clear

/*
  Construct table of descriptive statistics children born before the reform
  cut-off (by past CPS involvement)
*/
keep if post1 == 0
gen rev_cps_age16_binary = (cps_age16_binary == 0)

local var "enrol16 enrol17 cps_after16_binary evered16 evered_17to21 sex babywght adv_wgt babygest adv_b_gest adv_cong adv_ses mother_age adv_m_age adv_m_marstat"

eststo clear

eststo: estpost sum `var'
eststo: estpost sum `var' if cps_age16_binary == 0
eststo: estpost sum `var' if cps_age16_binary == 1
eststo: estpost ttest `var', by(rev_cps_age16_binary) unequal

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableA5.txt", replace label ///
	cells("mean(fmt(3) pattern (1 1 1 0)) b(star fmt(3) pattern (0 0 0 1))") ///
	mtitles("All" "No CPS" "Past CPS" "Delta, p < 0.01") collabels(none) ///
	float nogaps star(+ 0.01) varwidth(50) modelwidth(20) ///
	title("Appendix Table A5")
		
********************************************************************************
*	Appendix Table A6: Pre-reform census data on employment in South Australia
********************************************************************************
/*
  This table was constructed manually using aggregated data from the 2006 Census
  of Population and Housing. This aggregated data is available via the TableBuilder
  platform hosted by the Australian Bureau of Statistics (ABS) on the following
  website: https://www.abs.gov.au/statistics/microdata-tablebuilder/tablebuilder.
  To access the data, users will have to register a TableBuilder account with the
  ABS. Once registered, access the 2006 Census of Population of Housing dataset
  and select "Counting Persons, Place of Enumeration" as the unit of analysis.
  The relevant variables (employment status, hours worked, and weekly income) are
  available under the "Employment & Income" folder under "Person Variables".
  These can be filtered by education status and the statistics can then be exported
  into a .csv file. For hours worked and weekly income, we apply a mid-point method
  to the available categories in order to calculate the average hours worked and
  average weekly income. Note that the ABS applies perturbation to the data for 
  confidentiality purposes, so each user will generate a marginally different set
  of statistics.
*/

********************************************************************************
*	Appendix Table A7: RD estimates for leaving school at age 16 (by exit reason)
********************************************************************************

use "${clean}analysis_final", clear

*** Panel A: RD estimates for entire sample
est clear
local h = 721

foreach y in leave16_vet leave16_employment leave16_other{
	eststo: reghdfe `y' post1 h2 c.h2#post1 ///
		if rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth)
		
	qui sum `y' if rec_census15 == 1 & post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableA7.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("VET" "Employment" "Other") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast title("Appendix Table A7: Panel A")
	
*** Panel B: RD estimates by past CPS involvement
est clear

local h = 721

foreach y in leave16_vet leave16_employment leave16_other{
	eststo: reghdfe `y' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
		
	qui sum `y' if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Appendix_TableA7.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("VET" "Employment" "Other") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast title("Appendix Table A7: Panel B")

********************************************************************************
*	Appendix Table A8: RD estimates for new CPS notifications (public school sample)
********************************************************************************

use "${clean}analysis_final", clear

/*
  Restrict analysis to children in the public school sample (i.e. where we
  observe schooling outcomes)
*/
keep if rec_census15 == 1

* Construc table of RD estimates
*** Panel A: RD estimates for the entire sample
est clear
local h = 721

local vars cps_after16_binary cps_after16_phy cps_after16_sexual cps_after16_emo cps_after16_neg

foreach y in `vars'{
	qui eststo: reghdfe `y' post1 h2 c.h2#post1 ///
	if abs(h2) <= `h', vce(robust) absorb(monthofbirth)
	
	qui sum `y' if post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableA8.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("All" "Physical" "Sexual" "Emotional" "Neglect") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast title("Appendix Table A8: Panel A")

*** Panel B: RD estimates by past CPS involvement
eststo clear
local h = 721

local vars cps_after16_binary cps_after16_phy cps_after16_sexual cps_after16_emo cps_after16_neg

foreach y in `vars'{
	eststo: reghdfe `y' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Appendix_TableA8.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("All" "Physical" "Sexual" "Emotional" "Neglect") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast title("Appendix Table A8: Panel B")

********************************************************************************
*	Appendix Table B1: Robustness of RD estimates to alternative bandwidths
********************************************************************************

use "${clean}analysis_final", clear

*** Panel A: Optimal bandwidth (h = 721)
est clear
local h = 721

* Schooling outcomes
local varlist enrol16 enrol17

foreach z in `varlist'{
	eststo: reghdfe `z' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
}
		
* CPS outcome
eststo: reghdfe cps_after16_binary nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)

* ED outcomes
local varlist evered16 evered_17to21

foreach z in `varlist'{
	eststo: reghdfe `z' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
}

esttab using "${output}Appendix_TableB1.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) nostar ///
	float obslast modelwidth(25) ///
	title("Appendix Table B1: Panel A")

*** Panel B: One-and-a-half the optimal bandwidth (1.5h = 1081)
est clear
local h = 1081

* Schooling outcomes
local varlist enrol16 enrol17

foreach z in `varlist'{
	eststo: reghdfe `z' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
}
		
* CPS outcome
eststo: reghdfe cps_after16_binary nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)

* ED outcomes
local varlist evered16 evered_17to21

foreach z in `varlist'{
	eststo: reghdfe `z' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
}

esttab using "${output}Appendix_TableB1.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) nostar ///
	float obslast modelwidth(25) ///
	title("Appendix Table B1: Panel B")

*** Panel C: Half the optimal bandwidth (0.5h = 360)
/*
  Note that here, we are unable to include month-of-birth fixed effects directly
  in the regression model due to collinearity with the treatment/reform indicator
  (i.e. there are less than 12 months on either side of the cut-off). Therefore,
  we first run an auxiliary regression of the outcomes on month-of-birth across
  a 12-month period (~370 days)to account for seasonality. The residuals generated
  from this regression are then used as the outcome variable in the full RD 
  specification.
*/
est clear
local h = 360

* Schooling outcomes
local varlist enrol16 enrol17

foreach z in `varlist'{
	* Estimate auxiliary regression and generate residuals
	capture drop resid
	reg `z' i.monthofbirth#cps_age16_binary if rec_census15 == 1 & abs(h2) <= 370
	predict resid, residual
	
	* Use residuals as the outcome in the the RD specification
	local y resid

	eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h'  & rec_census15 == 1, robust
}

* CPS outcome
* Estimate auxiliary regression and generate residuals
capture drop resid
reg cps_after16_binary i.monthofbirth#cps_age16_binary if abs(h2) <= 370
predict resid, residual
	
* Use residuals as the outcome in the the RD specification
local y resid

eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
	c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
	if abs(h2) <= `h', robust

* ED outcomes
local varlist evered16 evered_17to21

foreach z in `varlist'{
	* Estimate auxiliary regression and generate residuals
	capture drop resid
	reg `z' i.monthofbirth#cps_age16_binary if abs(h2) <= 370
	predict resid, residual
	
	* Use residuals as the outcome in the the RD specification
	local y resid

	eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', robust
}

esttab using "${output}Appendix_TableB1.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) nostar ///
	float obslast modelwidth(25) ///
	title("Appendix Table B1: Panel C")
		
*** Panel D: Quarter the optimal bandwidth (0.25h = 180)
/*
  Note that here, we are unable to include month-of-birth fixed effects directly
  in the regression model due to collinearity with the treatment/reform indicator
  (i.e. there are less than 12 months on either side of the cut-off). Therefore,
  we first run an auxiliary regression of the outcomes on month-of-birth across
  a 12-month period (~370 days)to account for seasonality. The residuals generated
  from this regression are then used as the outcome variable in the full RD 
  specification.
*/
est clear
local h = 180

* Schooling outcomes
local varlist enrol16 enrol17

foreach z in `varlist'{
	* Estimate auxiliary regression and generate residuals
	capture drop resid
	reg `z' i.monthofbirth#cps_age16_binary if rec_census15 == 1 & abs(h2) <= 370
	predict resid, residual
	
	* Use residuals as the outcome in the the RD specification
	local y resid

	eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h'  & rec_census15 == 1, robust
}

* CPS outcome
* Estimate auxiliary regression and generate residuals
capture drop resid
reg cps_after16_binary i.monthofbirth#cps_age16_binary if abs(h2) <= 370
predict resid, residual
	
* Use residuals as the outcome in the the RD specification
local y resid

eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
	c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
	if abs(h2) <= `h', robust

* ED outcomes
local varlist evered16 evered_17to21

foreach z in `varlist'{
	* Estimate auxiliary regression and generate residuals
	capture drop resid
	reg `z' i.monthofbirth#cps_age16_binary if abs(h2) <= 370
	predict resid, residual
	
	* Use residuals as the outcome in the the RD specification
	local y resid

	eststo: reg `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', robust
}

esttab using "${output}Appendix_TableB1.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) nostar ///
	float obslast modelwidth(25) ///
	title("Appendix Table B1: Panel D")

********************************************************************************
*	Appendix Table B2: Robustness of RD estimates to using a quadratic control function
********************************************************************************

use "${clean}analysis_final.dta", clear

* Generate quadratic terms of the running variable
gen h2_squared = h2*h2
gen h3_squared = h3*h3

/*
  Perform optimal bandwidth selection for a quadratic RD model, using a cross-
  validation (CV) procedure which accounts for the quadratic running variable.
*/
* Define range of bandwidths under consideration (between 50 and 200 weeks)
local low 50
local high 200

* Generate counts for the corresponding row numbers (to store variables)
gen bandwidth = _n if inrange(_n, `low', `high')
	
* Set choice of window length for the CV procedure
local w 4 // 4 weeks or ~1 month
	
* Schooling outcomes
quietly{
foreach var in enrol16 enrol17{

	gen `var'_1 = 0

	forval bw = `low'/`high'{

	gen sqe = .

	forval x = 1/`w'{
		* Estimate regression model
		reg `var' nocps_post1 anycps_post1 cps_age16_binary h3 h3_squared ///
			c.h3#post1 c.h3_squared#post1 c.h3#cps_age16_binary c.h3_squared#cps_age16_binary ///
			c.h3#post1#cps_age16_binary c.h3_squared#post1#cps_age16_binary ///
			i.monthofbirth#cps_age16_binary ///
			if inrange(h3, -`bw' - `x', `bw' + `x') & !inrange(h3, -`x', `x') ///
			& rec_census15 == 1

		* Calculate prediction error in the `n' excluded months
		predict yhat_`var'
		replace sqe = (`var' - yhat_`var')^2 if (h3 == -`x' | h3 == `x') & rec_census15 == 1
		
		drop yhat_`var'
	}
	
	* Calculate mean-squared error (MSE) associated with the choice of bandwidth
	egen sse = sum(sqe)
	count if inrange(h3, -`w', `w') & h3 != 0 & rec_census15 == 1
	local N = r(N)
	replace `var'_1 = `var'_1 + (sse/(`N')) if bandwidth == `bw'
	
	drop sqe sse
}
}
}

* CPS outcome
quietly{
foreach var in cps_after16_binary{

	gen `var'_1 = 0

	forval bw = `low'/`high'{

	gen sqe = .

	forval x = 1/`w'{
		* Estimate regression model
		reg `var' nocps_post1 anycps_post1 cps_age16_binary h3 h3_squared ///
			c.h3#post1 c.h3_squared#post1 c.h3#cps_age16_binary c.h3_squared#cps_age16_binary ///
			c.h3#post1#cps_age16_binary c.h3_squared#post1#cps_age16_binary ///
			i.monthofbirth#cps_age16_binary ///
			if inrange(h3, -`bw' - `x', `bw' + `x') & !inrange(h3, -`x', `x')
		
		* Calculate prediction error in the `n' excluded months
		predict yhat_`var'
		replace sqe = (`var' - yhat_`var')^2 if (h3 == -`x' | h3 == `x')
		
		drop yhat_`var'
	}
	* Calculate mean squared error (MSE) associated with the choice of bandwidth
	egen sse = sum(sqe)
	count if inrange(h3, -`w', `w') & h3 != 0
	local N = r(N)
	replace `var'_1 = `var'_1 + (sse/(`N')) if bandwidth == `bw'
	
	drop sqe sse
}
}
}

* ED outcomes
quietly{
foreach var in evered16 evered_17to21{

	gen `var'_1 = 0

	forval bw = `low'/`high'{

	gen sqe = .

	forval x = 1/`w'{
		* Estimate regression model
		reg `var' nocps_post1 anycps_post1 cps_age16_binary h3 h3_squared ///
			c.h3#post1 c.h3_squared#post1 c.h3#cps_age16_binary c.h3_squared#cps_age16_binary ///
			c.h3#post1#cps_age16_binary c.h3_squared#post1#cps_age16_binary ///
			i.monthofbirth#cps_age16_binary ///
			if inrange(h3, -`bw' - `x', `bw' + `x') & !inrange(h3, -`x', `x')
			
		* Calculate prediction error in the `n' excluded months
		predict yhat_`var'
		replace sqe = (`var' - yhat_`var')^2 if (h3 == -`x' | h3 == `x')
		
		drop yhat_`var'
		}
		
	* Calculate mean-squared error (MSE) associated with the choice of bandwidth	
	egen sse = sum(sqe)
	count if inrange(h3, -`w', `w') & h3 != 0
	local N = r(N)
	replace `var'_1 = `var'_1 + (sse/(`N')) if bandwidth == `bw'
	
	drop sqe sse
	}
}
}

* Calculate the value of the joint CV criterion for each bandwidth
gen joint_cv = 0

local varlist enrol16 enrol17 cps_after16_binary evered16 evered_17to21

foreach var in `varlist'{
	replace joint_cv = joint_cv + `var'_1 
}

replace joint_cv = . if !inrange(bandwidth, 50, 200)

* Identify bandwidth which minimises the joint CV criterion
egen min_cv = min(joint_cv)
gen flag_min_bw = bandwidth if joint_cv == min_cv
egen min_bw = min(flag_min_bw)

* Now we use the optimal bandwidth (in days) to estimate the RD treatment effects
*** Panel A: RD estimates for entire sample
local h = min_bw*7
est clear

local vars enrol16 enrol17 cps_after16_binary evered16 evered_17to21

foreach y in `vars'{
    eststo: reghdfe `y' post1 h2 c.h2#post1 h2_squared c.h2_squared#post1 ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth)
		
	qui sum `y' if post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableB2.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(25) title("Appendix Table B2: Panel A")

*** Panel B: RD estimates by past CPS involvement
local h = min_bw*7
est clear

local vars enrol16 enrol17 cps_after16_binary evered16 evered_17to21

foreach y in `vars'{
    eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 h2_squared ///
		c.h2#post1 c.h2_squared#post1 c.h2#cps_age16_binary c.h2_squared#cps_age16_binary ///
		c.h2#post1#cps_age16_binary c.h2_squared#post1#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Appendix_TableB2.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(25) title("Appendix Table B2: Panel B")
	
********************************************************************************
*	Appendix Table B3: Robustness of RD estimates to using a global RD model
********************************************************************************

use "${clean}analysis_final.dta", clear

*** Panel A: RD estimates for entire sample
est clear

local vars enrol16 enrol17 cps_after16_binary evered16 evered_17to21

foreach y in `vars'{
    eststo: reghdfe `y' post1 h2 c.h2#post1, vce(robust) absorb(monthofbirth)
		
	qui sum `y' if post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableB3.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast varwidth(25) modelwidth(25) title("Appendix Table B3: Panel A")
	
*** Panel B: RD estimates by past CPS involvement
est clear

local vars enrol16 enrol17 cps_after16_binary evered16 evered_17to21

foreach y in `vars'{
    eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 ///
		c.h2#post1 c.h2#cps_age16_binary c.h2#post1#cps_age16_binary ///
		, vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Appendix_TableB3.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment at 16" "Enrolment at 17" "Any CPS notification" "Any ED visit at 16" "Any ED visit, 17 to 21") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(25) modelwidth(25) title("Appendix Table B3: Panel B")

********************************************************************************
*	Appendix Table B4: Robustness to alternative methods of inference
********************************************************************************

use "${clean}analysis_final", clear
program drop _all

* Generate quadratic terms of the running variable (for applying Noack and Rothe (2024))
gen h2_squared = h2*h2

/*
  Define program to calculate sharpened q-values following Anderson (2008). The
  sharpened q-values are calculated within domains of outcomes (schooling, 
  maltreatment, and ED visits) and within each subgroup (non-CPS vs any CPS).
  This program is intended to calculate the sharpened q-values for the non-CPS
  children.
*/
quietly{
program anderson2008_nocps

capture drop pval original_sorting_order rank bky06_qval
	* calculate sharpened q-value; taken from Michael Anderson's code (2008)

quietly gen float pval = .

display "***********************************"
display "Please paste the vector of p-values that you wish to test into the variable 'pval'"
display    "After pasting, type 'q' to resume"
display "***********************************"

forval j = 1/12{
    capture replace pval = pval_`j'_nocps if _n == `j'
}

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.

while `qval' > 0 {
    * First Stage
    * Generate the adjusted first stage q level we are testing: q' = q/1+q
    local qval_adj = `qval'/(1+`qval')
    * Generate value q'*r/M
    gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
    * Generate binary variable checking condition p(r) <= q'*r/M
    gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
    * Generate variable containing p-value ranks for all p-values that meet above condition
    gen reject_rank1 = reject_temp1*rank
    * Record the rank of the largest p-value that meets above condition
    egen total_rejected1 = max(reject_rank1)

    * Second Stage
    * Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
    local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
    * Generate value q_2st*r/M
    gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
    * Generate binary variable checking condition p(r) <= q_2st*r/M
    gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
    * Generate variable containing p-value ranks for all p-values that meet above condition
    gen reject_rank2 = reject_temp2*rank
    * Record the rank of the largest p-value that meets above condition
    egen total_rejected2 = max(reject_rank2)

    * A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
    replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
    * Reduce q by 0.001 and repeat loop
    drop fdr_temp* reject_temp* reject_rank* total_rejected*
    local qval = `qval' - .001
}
    
quietly sort original_sorting_order

end
}

/*
  Define program to calculate sharpened q-values following Anderson (2008). The
  sharpened q-values are calculated within domains of outcomes (schooling, 
  maltreatment, and ED visits) and within each subgroup (non-CPS vs any CPS).
  This program is intended to calculate the sharpened q-values for the children
  with past CPS involvement.
*/
quietly{
program anderson2008_anycps

capture drop pval original_sorting_order rank bky06_qval
	* calculate sharpened q-value; taken from Michael Anderson's code (2008)

quietly gen float pval = .

display "***********************************"
display "Please paste the vector of p-values that you wish to test into the variable 'pval'"
display    "After pasting, type 'q' to resume"
display "***********************************"

forval j = 1/12{
    capture replace pval = pval_`j'_anycps if _n == `j'
}

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.

while `qval' > 0 {
    * First Stage
    * Generate the adjusted first stage q level we are testing: q' = q/1+q
    local qval_adj = `qval'/(1+`qval')
    * Generate value q'*r/M
    gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
    * Generate binary variable checking condition p(r) <= q'*r/M
    gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
    * Generate variable containing p-value ranks for all p-values that meet above condition
    gen reject_rank1 = reject_temp1*rank
    * Record the rank of the largest p-value that meets above condition
    egen total_rejected1 = max(reject_rank1)

    * Second Stage
    * Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
    local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
    * Generate value q_2st*r/M
    gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
    * Generate binary variable checking condition p(r) <= q_2st*r/M
    gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
    * Generate variable containing p-value ranks for all p-values that meet above condition
    gen reject_rank2 = reject_temp2*rank
    * Record the rank of the largest p-value that meets above condition
    egen total_rejected2 = max(reject_rank2)

    * A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
    replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
    * Reduce q by 0.001 and repeat loop
    drop fdr_temp* reject_temp* reject_rank* total_rejected*
    local qval = `qval' - .001
}
    
quietly sort original_sorting_order

end
}

/*
  Now we construct the table of clustered standard errors, honest p-values (using
  the 'rdhonest' package in Stata), and sharpened q-values (using the program
  defined above based on code by Michael Anderson). Note that the 'rdhonest'
  package's functionality for including covariates does not work properly, and
  so we so adjust for the month-of-birth fixed effects manually by applying the
  Frisch-Waugh-Lovell (FWL) theorem.
*/
***** Schooling outcomes: Columns (1) and (2)
est clear

local h = 721
local x = 0
local k = 0

foreach y in enrol16 enrol17{
    local x = `x' + 1
	local k = `k' + 1
	
	*** Clustered standard errors
	eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if rec_census15 == 1 & abs(h2) <= `h', vce(cluster h2) absorb(monthofbirth#cps_age16_binary)
		
	mat m2 = r(table)
	scalar cluster_nocp_p = m2[4,1]
	scalar cluster_anycp_p = m2[4,2]
	
	mat cluster_p = (cluster_nocp_p, cluster_anycp_p)
	mat colnames cluster_p = nocp_post1 anycp_post1
	estadd matrix cluster_p: est`x'
	
	*** Honest p-values
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the non-CPS children (i.e. Reform x No CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 0 & rec_census15 == 1, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 0 & rec_census15 == 1 ///
		& abs(h2) <= 721
	gen adj_nocp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_nocp_`y' h2 if cps_age16_binary == 0 & rec_census15 == 1, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_nocp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the past CPS children (i.e. Reform x Past CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 1 & rec_census15 == 1, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 1 & rec_census15 == 1 ///
		& abs(h2) <= 721
	gen adj_anycp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_anycp_`y' h2 if cps_age16_binary == 1 & rec_census15 == 1, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_anycp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	* Store rdhonest p-values in matrix (to export later)
	mat rdhonest_p = (rdhonest_nocp_p, rdhonest_anycp_p)
	mat colnames rdhonest_p = nocps_post1 anycps_post1
	estadd matrix rdhonest_p: est`x'
	
	*** Store p-values for calculating sharpened q-values
	reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if  rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	mat m2 = r(table)
	scalar pval_`k'_nocps = m2[4,1]
	scalar pval_`k'_anycps = m2[4,2]
}

/*
  Now we use the earlier defined programs to calculate sharpened q-values,
  separately for the non-CPS children and the past CPS children.
*/
anderson2008_nocps
	
forval z = 1/2{
	scalar qval_`z'_nocps = bky06_qval[`z']
}

anderson2008_anycps

forval z = 1/2{
	scalar qval_`z'_anycps = bky06_qval[`z']
}

* Store sharpened q-values in matrix (to export later)
local x = 0
forval z = 1/2{
    local x = `x' + 1
	
    mat qval_p = (qval_`z'_nocps, qval_`z'_anycps)
	mat colnames qval_p = nocps_post1 anycps_post1
	estadd matrix qval_p: est`x'
}

***** CPS outcomes: Columns (3) to (7)
local h = 721
local vars cps_after16_binary cps_after16_phy cps_after16_sexual cps_after16_emo cps_after16_neg
scalar drop _all

local x = 2
local k = 0

foreach y in `vars'{
    local x = `x' + 1
	local k = `k' + 1
	
	*** Clustered standard errors
	eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(cluster h2) absorb(monthofbirth#cps_age16_binary)
		
	mat m2 = r(table)
	scalar cluster_nocp_p = m2[4,1]
	scalar cluster_anycp_p = m2[4,2]
	
	mat cluster_p = (cluster_nocp_p, cluster_anycp_p)
	mat colnames cluster_p = nocps_post1 anycps_post1
	estadd matrix cluster_p: est`x'
	
	*** Honest p-values
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the non-CPS children (i.e. Reform x No CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 0, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 0 & abs(h2) <= 721
	gen adj_nocp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_nocp_`y' h2 if cps_age16_binary == 0, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_nocp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))

	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the past CPS children (i.e. Reform x Past CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 1, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 1 & abs(h2) <= 721
	gen adj_anycp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_anycp_`y' h2 if cps_age16_binary == 1, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_anycp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	* Store honest p-values in matrix (to export later)
	mat rdhonest_p = (rdhonest_nocp_p, rdhonest_anycp_p)
	mat colnames rdhonest_p = nocps_post1 anycps_post1
	estadd matrix rdhonest_p: est`x'
	
	*** Store p-values for calculating sharpened q-values
	reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	mat m2 = r(table)
	scalar pval_`k'_nocps = m2[4,1]
	scalar pval_`k'_anycps = m2[4,2]
}

/*
  Now we use the earlier defined programs to calculate sharpened q-values,
  separately for the non-CPS children and the past CPS children.
*/
anderson2008_nocps
	
forval z = 1/5{
	scalar qval_`z'_nocps = bky06_qval[`z']
}

anderson2008_anycps

forval z = 1/5{
	scalar qval_`z'_anycps = bky06_qval[`z']
}

* Store sharpened q-values in matrix (to export later)
local x = 2
forval z = 1/5{
    local x = `x' + 1
	
    mat qval_p = (qval_`z'_nocps, qval_`z'_anycps)
	mat colnames qval_p = nocps_post1 anycps_post1
	estadd matrix qval_p: est`x'
}

***** All-cause ED outcomes: Columns (8) and (9)
local h = 721
local vars evered16 evered_17to21
scalar drop _all

local x = 7
local k = 0

foreach y in `vars'{
    local x = `x' + 1
	local k = `k' + 1
	
	*** Clustered standard errors
	eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(cluster h2) absorb(monthofbirth#cps_age16_binary)
		
	mat m2 = r(table)
	scalar cluster_nocp_p = m2[4,1]
	scalar cluster_anycp_p = m2[4,2]
	
	mat cluster_p = (cluster_nocp_p, cluster_anycp_p)
	mat colnames cluster_p = nocps_post1 anycps_post1
	estadd matrix cluster_p: est`x'
	
	*** Honest p-values
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the non-CPS children (i.e. Reform x No CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 0, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 0 & abs(h2) <= 721
	gen adj_nocp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_nocp_`y' h2 if cps_age16_binary == 0, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_nocp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the past CPS children (i.e. Reform x Past CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 1, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 1 & abs(h2) <= 721
	gen adj_anycp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_anycp_`y' h2 if cps_age16_binary == 1, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_anycp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	* Store honest p-values in matrix (to export later)
	mat rdhonest_p = (rdhonest_nocp_p, rdhonest_anycp_p)
	mat colnames rdhonest_p = nocps_post1 anycps_post1
	estadd matrix rdhonest_p: est`x'
	
	*** Store p-values for calulating sharpened q-values
	reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	mat m2 = r(table)
	scalar pval_`k'_nocps = m2[4,1]
	scalar pval_`k'_anycps = m2[4,2]
}

/*
  Now we use the earlier defined programs to calculate sharpened q-values,
  separately for the non-CPS children and the past CPS children.
*/
anderson2008_nocps
	
forval z = 1/2{
	scalar qval_`z'_nocps = bky06_qval[`z']
}

anderson2008_anycps

forval z = 1/2{
	scalar qval_`z'_anycps = bky06_qval[`z']
}

* Store sharpened q-values in matrix (to export later)
local x = 7
forval z = 1/2{
    local x = `x' + 1
	
    mat qval_p = (qval_`z'_nocps, qval_`z'_anycps)
	mat colnames qval_p = nocps_post1 anycps_post1
	estadd matrix qval_p: est`x'
}

***** Cause-specific ED outcomes: Columns (10) to (13)
local h = 721
local vars everedinj16 everedment16 evereddigest16 everedoth16
scalar drop _all

local x = 9
local k = 0

foreach y in `vars'{
    local x = `x' + 1
	local k = `k' + 1
	
	*** Clustered standard errors
	eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(cluster h2) absorb(monthofbirth#cps_age16_binary)
		
	mat m2 = r(table)
	scalar cluster_nocp_p = m2[4,1]
	scalar cluster_anycp_p = m2[4,2]
	
	mat cluster_p = (cluster_nocp_p, cluster_anycp_p)
	mat colnames cluster_p = nocps_post1 anycps_post1
	estadd matrix cluster_p: est`x'
	
	*** Honest p-values
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the non-CPS children (i.e. Reform x No CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 0, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 0 & abs(h2) <= 721
	gen adj_nocp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_nocp_`y' h2 if cps_age16_binary == 0, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_nocp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	/*
	  Calculate rule-of-thumb for 'm' following Noack and Rothe (2024), apply FWL
	  to adjust for month-of-birth fixed effects, and then calculate honest p-values
	  for the past CPS children (i.e. Reform x Past CPS) using the 'rdhonest' package.
	*/
	reg `y' h2 h2_squared if post1 == 0 & cps_age16_binary == 1, robust
	local m = abs(_b[h2_squared]*2)
	
	reg `y' post1 h2 c.h2#post1 mb* if cps_age16_binary == 1 & abs(h2) <= 721
	gen adj_anycp_`y' = `y' - _b[mb1]*mb1 - _b[mb2]*mb2 - _b[mb3]*mb3 - _b[mb4]*mb4 ///
		- _b[mb5]*mb5 - _b[mb6]*mb6 - _b[mb7]*mb7 - _b[mb8]*mb8 - _b[mb9]*mb9 ///
		- _b[mb10]*mb10 - _b[mb11]*mb11
	
	rdhonest adj_anycp_`y' h2 if cps_age16_binary == 1, m(`m') kernel("uni") h(721)
	
	local b = e(bias)/e(se)
	local t = e(est)/e(se)
	scalar rdhonest_anycp_p = normal(`b' - abs(`t')) + normal(-`b' - abs(`t'))
	
	* Store honest p-values in matrix (to export later)
	mat rdhonest_p = (rdhonest_nocp_p, rdhonest_anycp_p)
	mat colnames rdhonest_p = nocps_post1 anycps_post1
	estadd matrix rdhonest_p: est`x'
	
	*** Store p-values for calculating sharpened q-values
	reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	mat m2 = r(table)
	scalar pval_`k'_nocps = m2[4,1]
	scalar pval_`k'_anycps = m2[4,2]
}

/*
  Now we use the earlier defined programs to calculate sharpened q-values,
  separately for the non-CPS children and the past CPS children.
*/
anderson2008_nocps
	
forval z = 1/4{
	scalar qval_`z'_nocps = bky06_qval[`z']
}

anderson2008_anycps

forval z = 1/4{
	scalar qval_`z'_anycps = bky06_qval[`z']
}

* Store sharpened q-values in matrix (to export later)
local x = 9
forval z = 1/4{
    local x = `x' + 1
	
    mat qval_p = (qval_`z'_nocps, qval_`z'_anycps)
	mat colnames qval_p = nocps_post1 anycps_post1
	estadd matrix qval_p: est`x'
}

/*
  Now we can export the entire table using the 'esttab' package
*/
esttab using "${output}Appendix_TableB4.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps ///
	mtitles("Enrolment: 16" "Enrolment: 17" "CPS: all" "CPS: physical" "CPS: sexual" "CPS: emotional" "CPS: neglect" "Any ED: 16" "Any ED: 17 to 21" "ED: injury" "ED: mental" "ED: digestive" "ED: other") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	float nostar collabels(none) ///
	cells(b(fmt(%9.3f)) se(fmt(%9.3f) par pvalue(cluster_p)) ///
	rdhonest_p(fmt(%9.3f) par([ ]) pvalue(rdhonest_p)) ///
	qval_p(fmt(%9.3f) par({ })  pvalue(qval_p))) ///
	modelwidth(10) ///
	title ("Appendix Table B4")
	
********************************************************************************
*	Appendix Table B5: Placebo analysis
********************************************************************************

use "${clean}analysis_final.dta", clear

/*
  Generate running variable for placebo analysis (i.e. day-of-birth from the 
  placebo cut-off, 1 January 1991)
*/
gen h2_placebo = dob - mdy(1, 1, 1991)

* Generate placebo reform/treatment indicators (i.e. born after 1 January 1991)
gen post1_placebo = (h2_placebo >= 0)
gen nocps_post1_placebo = (h2_placebo >= 0 & cps_age16_binary == 0)
gen anycps_post1_placebo = (h2_placebo >= 0 & cps_age16_binary == 1)

* Generate RD estimates using the placebo cut-off
est clear

local h = 650

local vars cps_after16_binary evered16

foreach y in `vars'{
	eststo: reghdfe `y' nocps_post1_placebo anycps_post1_placebo cps_age16_binary ///
	h2_placebo c.h2_placebo#post1_placebo ///
	c.h2_placebo#cps_age16_binary c.h2_placebo#post1_placebo#cps_age16_binary ///
	if abs(h2_placebo) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	qui sum `y' if post1_placebo == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if post1_placebo == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableB5.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("Any CPS notification" "Any ED visit at age 16") ///
	keep(nocps_post1_placebo anycps_post1_placebo) ///
	coeflabel(nocps_post1_placebo "Reform x No CPS" ///
	anycps_post1_placebo "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Placebo pre-reform mean (no CPS)" ///
	"mean_anycp Placebo pre-reform mean (past CPS)") ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	float obslast varwidth(40) modelwidth(30) title("Appendix Table B5")
	
********************************************************************************
*	Appendix Table B6: RD model including additional adversity indicators
********************************************************************************

use "${clean}analysis_final.dta", clear

/*
  Generate treatment/reform indicators interacted with mother's age, mother's 
  marital status (married/de facto), and socioeconomic status (i.e. low SES) at birth
*/
foreach x in adv_m_age adv_m_marstat adv_ses{
	gen no`x'_post1 = 1 if h2 >= 0 & `x' == 0
	replace no`x'_post1 = 0 if h2 < 0 | `x' == 1

	gen any`x'_post1 = 1 if h2 >= 0 & `x' == 1
	replace any`x'_post1 = 0 if h2 < 0 | `x' == 0
}
	
* Construct table of RD estimates (using RD model with new interaction terms)
eststo clear

local h = 721
local vars cps_after16_binary evered16

foreach y in `vars'{
	eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post1 ///
		c.h2#cps_age16_binary c.h2#post1#cps_age16_binary ///
		anyadv_ses_post1 adv_ses c.h2#adv_ses c.h2#post1#adv_ses ///
		anyadv_m_marstat_post1 adv_m_marstat c.h2#adv_m_marstat c.h2#post1#adv_m_marstat ///
		anyadv_m_age_post1 adv_m_age c.h2#adv_m_age c.h2#post1#adv_m_age ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
			
	qui sum `y' if post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
		
	qui sum `y' if post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableB6.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("Any CPS notification" "Any ED visit") ///
	keep(nocps_post1 anycps_post1 anyadv_ses_post1 anyadv_m_marstat_post1 anyadv_m_age_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS" ///
	anyadv_ses_post1 "Reform x Low SES" ///
	anyadv_m_marstat_post1 "Reform x Mother not married/de facto" ///
	anyadv_m_age_post1 "Reform x Mother's age < 21") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///
	star(* 0.05 ** 0.01) ///
	float obslast varwidth(40) modelwidth(25) title("Appendix Table B6")
	
********************************************************************************
*	Appendix Table C1: RD estimates of the reform's effect on school absenteeism
********************************************************************************

use "${clean}analysis_final.dta", clear

* Construct table of RD estimates
est clear

local h = 721

eststo: reghdfe days_present16 nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
	c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
	if rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
qui sum days_present16 if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 0
estadd scalar mean_nocp = r(mean)
	
qui sum days_present16 if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 1
estadd scalar mean_anycp = r(mean)

* Export table using the 'esttab' package
esttab using "${output}Appendix_TableC1.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("Days of school attended at age 16") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(40) title("Appendix Table C1")

/********************************************************************************
References:

Anderson, Michael L., "Multiple Inference and Gender Differences in the Effects of
	Early Intervention: A Reevaluation of the Abecedarian, Perry Preschool, and Early
	Training Projects," Journal of the American Statistical Association, 2008, 103 (484),
	1481–1495.
Gnanamanickam, Emmanuel S., Ha Nguyen, Jason M. Armfield, James C. Doidge, Derek S. Brown,
	David B. Preen, and Leonie Segal, "Child Maltreatment and Emergency
	Department Visits: A Longitudinal Birth Cohort Study from
	Infancy to Early Adulthood," Child Abuse & Neglect, 2022, 123, 105397.
Noack, Claudia and Christoph Rothe, "Bias-Aware Inference in Fuzzy Regression
	Discontinuity Designs," Econometrica, 92: 687-711.
********************************************************************************/